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Abstract. We study the influence of network topology on retrieval properties of 
recurrent neural networks, using replica techniques for diluted systems. The theory 
is presented for a network with an arbitrary degree distribution p{k) and applied to 
power law distributions p{k) ^ k~'^, i.e. to neural networks on scale-free graphs. A 
bifurcation analysis identifies phase boundaries between the paramagnetic phase and 
either a retrieval phase or a spin glass phase. Using a population dynamics algorithm, 
the retrieval overlap and spin glass order parameters may be calculated throughout the 
phase diagram. It is shown that there is an enhancement of the retrieval properties 
compared with a Poissonian random graph. We compare our findings with simulations. 
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1. Introduction 

The impressive ability of human and animal brains to recognize and manipulate complex 
patterns under real- world (i.e. noisy and often conflicting) conditions continues to appeal 
not only to biologists but also to physicists, computer scientists and engineers, albeit 
with the latter driven by different objectives and motivations. Hopfield [Tj was one of 
the first to introduce a simple model to describe associative memory in recurrent neural 
networks successfully, based on the biologically motivated Hebbian rule for adapting 
the connections between the neurons (the 'synapses'). His model initiated a period of 
intense research activity. The success of these early neural network models was mainly 
due to their analytic tractability, which was achieved upon sacrificing biological realism; 
all neurological connectivity structures were sacrificed by the first generation of (fully 
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connected) models. However, the research area was thereby able to benefit significantly 
from recent advanced in mean-field spin glass theory IllEJi and many new results were 
pubhshed in the second half of the eighties; see e.g. [UEIIE] or [Zj. 

A step towards increased biological realism was made by the introduction of diluted 
neural network models. Initially, in the thermodynamic limit each neuron was on average 
connected to a vanishing fraction of the system, but this fraction contained an infinite 
number of nodes. These models were solvable by virtue of the specific nature of their 
architectures: one either chooses strictly symmetric dilution (so detailed balance and 
hence equilibrium analysis are preserved, e.g. [311011111), or strictly asymmetric dilution, 
which ensures that neuron states are statistically independent on finite times jSl (now 
the local fields are described by Gaussian distributions, leading to simple dynamic order 
parameter equations). In the early models, the bond statistics were uniform over the 
entire network, leading to thin tails in its degree distribution, whereas the connectivity of 
a real neuron is known to vary strongly within the brain ^2] • In response to this, there 
have been several recent studies of recurrent neural network models with alternative 
connectivity distributions. Most evolve around numerical simulations of Hopfield-type 
models on graphs with power law degree distributions (THl (HI 113 HE]- Examples of 
recent analytic work on recurrent neural networks with finite connectivity can be found 
in [T71 UHl; both deal with Poissonian graphs and apply the equihbrium statistical 
mechanical techniques of dilute disordered spin systems [ini EUl 121] • 

The objective of this paper is to extend and generalize the solution for finitely 
connected Poissonnian neural networks ^3 ^| to recurrent neural networks with 
arbitrary degree distribution, in the spirit of |221 1211 within the replica-symmetric 
(RS) ansatz. We derive analytically phase diagrams for networks with Hebbian synapses 
and arbitrary degree distributions p{k), and apply population dynamics algorithms to 
obtain the values of the order parameters in the three phases (viz. paramagnetic, 
retrieval, and spin-glass). This study thereby establishes a connection between the 
equilibrium statistical mechanics of neural networks and the theory of so-called 'complex 
networks'. In line with biological reality, we find that recurrent neural networks with 
degree distributions with 'fat tails' are indeed superior to those with Poissonnian degree 
distributions, in terms of the size of the recall region in the phase diagram. 

2. Model definitions 

Our model is a system of Ising spin neurons G {—1,1}, with i = 1,...,N. 
The neurons are located on the nodes of a graph with arbitrary degree distribution 
p{k) = N~'^ J2i ^k,ki, where fcj denotes the number of neurons connected to neuron i. This 
system is assumed to be in thermodynamic equilibrium, described by the Hamiltonian 

N p N 

n = - Y: ^^J^J^-T.h''T.^^^^ (l) 

i<j=l M=l 

Here the {h^} represent generating fields, and the vectors (^f, . . . ,^^) represent stored 
random iV-bit patterns. We will abbreviate the bits to be stored at a given node i as 
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~ i^iy • • • 5^f)- Since the number of connections per neuron is finite, tlie number of 
patterns p must be of order 0{N^). The bonds Jjj depend on the patterns via 



J i' 



Cij 
Jk) 



(2) 



X, 



with ■ = Ylli=i^ii'j and {k) = J2k>(}P{k)k- Special cases of interest are (f){x) 
viz. Hebbian bonds, and (f){x) = sign(a;), viz. chpped Hebbian bonds. 

The variables Cij G {0, 1} specify our graph microscopically; we extend their 



definition to all pairs {i,j) by putting q 



and Cii = 0. It is known from complex 



network theory [211 ESI ^^^^ ^ connectivity distribution p{k) alone does not fully specify 
the statistics of a graph. Here we draw the matrix c = {cjj}, which represents quenched 
disorder for the spin system (P), randomly from the probability distribution 



V{c) 



with the single-bond probabilities 







Yii ^ki,^^^^ Cij 











Pic. 



Or 1 



Disorder averages ((^(c)))^ over the ensemble of graphs are thus given by |221I2S] 



c i<j i •' 



^ =E[nnQ,)^.,j[n/ 

c i<j i 



27r 



(3) 

(4) 

(5) 
(6) 



Another statistical quantity to characterize graphs, beyond is the degree-degree 

correlation: the joint probability u{ki,kj) that a pair of nodes i and j are connected, 
and have connectivities ki and kj, respectively. For the present ensemble (jS)) one finds 

p{ki)p{kj)kikj 



(jJ ^ ki "I kj 



{k)N 



(7) 



3. Replica calculation of the free energy and order parameters 

We calculate the free energy per spin and the relevant order parameters, using the replica 
techniques as developed for constrained connectivity graphs, along the lines of [21 
Thus the asymptotic free energy per spin / = — lim7v^oo(/3iV)~^ ^ogZ is written as 



/ = - hm hm log((Z")), 

N^oon^O J\lnp 

where -Z" is the usual n-replicated partition function 



N n 



^'^ = E---Ee^p [/?EE^r^^K + /5E^"EEer^r 

<T 1 <T jv a=l i<j fi=l 



i=l a=l 



(8) 



(9) 
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and cTj = [a], . . . , a") is the n-replicated spin at site i. Upon performing the trace over 
the Cij (i.e. the disorder average) one obtains 



Af 



N 

n 

i=l 



X exp 



E 

{k)_ 

2N , 



271 



-iipiki 



N n 



exp [pY.h'-EE^'^r 



fj,=l i=l a=l 



N 

J2 {e^^' 



Tfc) vvl,-^,)o-,-<T, _i\ ^ 0(iV°)l (10) 



J=l 



In order to go to an effective single-site problem one first introduces the concept of 
sublattices = {i\$,i = ^} and defines the following order parameter functions: 



These are reminiscent of the replicated spin probability distributions within sublattices, 
as in ^7] , but here include extra phase factors e*'^' whose effect is to replace ki by fcj — 1 
in expressions of the type (0). At the physical saddle-point in the subsequent calculation 
one finds the physical meaning of (lllj) to be 



(12) 



with (...) denoting a thermal average over the n-replicated spins. Thus the order 
parameter (fTTj) is just the distribution of a replicated cavity spin in sublattice ^. 
After some straightforward manipulations (factorization over sites, integration over the 
variables ipi, etc) we find / = lim„^o extr|pp|/[{P, P}], where 



{k) 



E((P^(^)P^'(-')e^^^^-^'^^"-"'^))^,^' 



(k) 



(13) 



2Pn^r ' ' "^'^ 2^n 

This involves the sublattice averages (/(C))^ = Z]^p(C)/(C)) withp^ = limAr_^oo l-^^l/A^ 
Varying (fT^ with respect to {P^(cr), P^(cr)} leads to the saddle point equations 



a' 



p,H = E ^ 



(k) e'^h-^I^^=i""P|(cr) 
Upon adopting the usual rephca symmetric (RS) ansatze [21111211122], i-e. 



(14) 
(15) 



P^i^) 



dh W^{h) 
du 



[2 cosh{(3h)f 
[2 cosh(/?M)]" 



(16) 
(17) 
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one finds the saddle-point equations (fT^ and (fT^ reducing to 

Q^(u) = {Jdh W^ih) 5 u- ^tanh^i [tanh(/3/i) tanh(^^^|--^)] 



k-l 



fc-i 



(18) 



(19) 



where W^{h) and Q^{u) are probabihty distributions for the effective cavity fields and 
propagated fields (messages) in sublattice ^, respectively ^23,. Insertion of (|16I17|) into 
(fT^ gives us the RS free energy per spin: 



Hi- 



(3f = {k){Jdhdu W^{h)Q ^{u) log [l + tanh{f3u)ta.nh{f3h)])^ 

-^{{J dhdh'W^{h)W^'{h') log [l + tanhiPh) tanh(/?/i') tanh [ ^'^^^^j ^ ^ 



J2Pk{ [Y[du£Q^{ui)^ log 



nii2cosh(/3u,) 



■((log cosh [ — - {k) log2 



(20) 



2 (fc) 

We may finally use the generating fields /i^ to find explicit expressions for the disorder- 
averaged pattern recall overlaps = limAr^oo Z^i C^(((<^i)))c = ''{df /dh'^)\h=o- 

k k 



Y.Pk f[\{du^QAu,)\ tanh(/3^w,)). 

k ■' £=1 e=i 



(21) 



In a similar manner one may derive an expression for the disorder-averaged RS spin- 
glass order parameter q = limN^oo J2i ((('^i)^))c; upon adding a term of the form 
XJ2a<p(^aC/3 to the replicated Hamiltonian in Q. The result is 

k k 

(22) 



(l = Y.Pk { [lidueQ^iui)^ tanh^ 



These expressions have a transparent interpretation, given that within the cavity 
formalism the local magnetization at a site i is indeed given by = tanh {f3 ue) ■ 



4. Phase diagram and order parameters 

The paramagnetic (P) phase, where (pT|l and (j22I) are zero, has W^{x) = Q^{x) = 6{x) 
for all ^. The recall phase (R) is defined by m'^ ^ for some /i. In the spin-glass (SG) 
phase, g > but m'^ = for all ii. Transitions away from the P phase are expected to 
be second order, allowing us to find the P— >(R,SG) transitions via a simple continuous 
bifurcation analysis. In contrast, locating the SG— >-R transition requires knowledge 
of the (nontrivial) functions W^{h) and Q^{u) in the R or SG regimes; to find these 
transitions we will solve p8ll9p numerically using a population dynamics algorithm. 

Following [13 HE] we apply a bifurcation analysis to compute the second order 
transition lines away from the paramagnetic phase. By assuming the effective fields to 
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be small, i.e. Jdh W^{h)h^ = O(e^) for all and expanding equation (|18I19|) up to 
order e^, one finds the following bifurcation conditions for transitions away from P: 



P^R: dhW^{h)h 



{k) 



i^(tanh[-^^^^-^') 



P^SG: dhW^{h)h 



(k) 



(k) 

(tanh-^^(^-^') 



{k) 



dhW^{h)h)^ (23) 



dh W^{h)h^)^ (24) 



For = 2 P (random patterns) one can find the eigenvalues of the relevant matrices 



M, 



Q 



jk') - (k) 

(k) 
jk') - {k) 

(k) 



tanh 



tanh 



ik) 



ik) 



and the conditions for a second order transition become, see jT7| : 



R : 



SG 



{k') - (k) 2 




tanh 



2n) tanh 



pcpip - 2n) 



(k) 



2n) 



(k) 



(25) 
(26) 

(27) 
(28) 



For a Poissonian degree distribution one has {k'^) = + (fc), and we recover the results 
in ^Tj. For non-Poissonian distributions, equation ()27|) predicts an enlargement of the 
retrieval phase when the degree distribution has fat tails, e.g. for power-law distributions 
p{k) ~ k~'^ . As noted by |2n|, if the second moment of the degree distribution is not 
finite (e.g. for power- law distributions with 7 < 3), there is always a retrieval phase in 
the thermodynamic limit (unless one re-scales temperature first). For a single pattern, 
the retrieval phase boundary is given, in accordance with by 



m 



(29) 



(our expression differs slightly, due to our rescaling of the bonds by a factor {k)~ 



To find our order parameters we solve equations ()18pi9|l numerically, using a 
population dynamics algorithm [211 I2H]- We have two distributions for each of our 
2^ sublattices, so the required CPU time grows exponentially with p. However, one 
may exploit sublattice symmetries, especially when the system is condensed in a single 
pattern retrieval state 29 . If the first pattern is condensed, numerical and analytical 
evidence (as yet short of a proof) suggests the solution of ()18pi9j) to have W^{h) = 
W^i{h). Moreover, for random patterns one expects the symmetry W^i{h) = W{C,^h). 
Insertion of this ansatz (and a similar one for Q^) into (jl8ll9|) leads us to 

fc-l k-l 



W{h) 
Q{u) 



kpk 

rik)j 



2P 



X 



[l[du, Q{u,)]6[h 
''dh W{h) 



£=1 



(30) 



n=0 




^[m — — tanh ^ [tanh(/3/i) tanh(^^^^— — — 



(k) 



(31) 
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Figure 1. RS phase diagrams for p{k) ~ k^'' with p{k) — ior k < m, where 7 — 3.1 
(left) and "f — 4 (right). Phase boundaries found by bifurcation analysis, separating 
paramagnetic (P) from retrieval (R) and spin glass (SG) phases. For 7 = 4 the R^SG 
phase boundary is calculated for T values which are multiples of 0.2, and is found via 
population dynamics (with a condensed ansatz). For 7 = 3.1, average connectivity 
values corresponding to m = 1 (circles), m = 2 (squares) and m — 3 (diamonds) are 
(fc) = 1.318, (k) = 3.055 and (k) — 4.898 respectively. For 7 = 4, average connectivity 
values corresponding to m = 2 (circles), m = 3 (squares) and m = 4 (diamonds) are 
(k) = 2.454, (fc) = 3.887 and (k) = 5.352 respectively. 



Inserting (j!?T|l into ()3()|1 then gives the relatively simple expression 



W{h) 



p-i 




ni=0 



1 



X 5<h — —y^^ tanh ^ tanh(/3/i£) tanh[ 



{k) 



(32) 



5. Comparison with simulations 



In figure ^ we present the resulting RS phase diagrams in the {a, T)-plane, where 
a = p/ {k) for Hopfield-type networks with power-law degree distribution p{k) ~ k~'^. 
The distributions p{k) are characterized by 7 = 3.1 and 7 = 4 respectively, and by a 
variable m which defines a lower cutoff {p{k) = for k < m). At high T one finds the 
paramagnetic (P) phase. At sufficiently low T one finds a retrieval (R) phase (small a) 
or a spin glass (SG) phase (large a). The value 7 = 3.1 is close to the critical value 
7c = 3 below which there is no paramagnetic phase, yet here the phase diagram is found 
to be similar to that corresponding to Poissonian graphs ^Tj (with re-scaled values of 
T and a). In the 7 = 4 phase diagram we also indicate the location of the R— >-SG 
transition, resulting from a population dynamics calculation. The P— >^R and P— s>SG 
boundaries were found by solving the bifurcation equations ()27I28|) numerically. For 



small values of the order parameters, i.e. large T or a, finding accurate numerical values 




Figure 2. Condensed recall overlap m as calculated from our RS theory via 
population dynamics (squares) versus numerical simulations (circles), for a graph with 
p{k) [k{k + l)(fc + 2)]"i (left), and a graph with p{k) - k^^ (right), both with 
p{k) — for fc < 3. The corresponding average connectivities are (fc) = 6 (left) and 
(k) — 5.123 (right). All data shown are averages over 10 runs, in both simulations 
(iV = 10'' spins) and in population dynamics (with populations of size N = 10^). 



for the order parameters becomes increasingly difficult. Hence, for large fluctuations in 
the connectivity (i.e. large values of (/c^)), one cannot expect accurate results for the 
R^SG transition on the basis of a population dynamics algorithm. This is why we have 
omitted the R— >SG line in the 7 = 3.1 phase diagram. 

In figure El we show as function of temperature the recall overlap (for states 
where only pattern 1 is condensed) as obtained from a population dynamics calculation, 
together with the measurements of in numerical simulations. Both simulation and 
population dynamics had N = lO'^ and all data are averages over 10 runs. The left figure 
refers to p{k) ~ [k{k + l)(/c + 2)]~^, corresponding to the degree distribution resulting 
from a Barabasi-Albert algorithm for network growth [SHI 1^- It should be noted, 
however, that the degree-degree correlation generated by the latter algorithm differs 
from the one ((Zj) in our present model (see e.g. [24,j)- Consequently, a different algorithm, 
similar to the one described in had to be used here. First, a set {ki, . . . , k^} is 
generated in accordance with the distribution p{k). One then chooses randomly two 
sites i and j with probabilities p(i) = ~ p{ki)ki. Unless these two sites coincide 

or already share a bond, they are connected. If site i already has ki connections, it is 
excluded from the process to speed up the algorithm. This process is repeated until all 
connectivities have the correct value. We also used this algorithm to generate power law 
distributed graphs with 7 = 3 (right figure). For both architectures, the overlaps m} 
have been plotted for networks with p = 1, p = 3 and p = 5. At low temperature, the 
results of the simulation agree well with the population dynamics results. For the values 
of p used in figure (j2I) one has pattern retrieval at sufficiently low temperatures. In fact 
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Figure 3. Average retrieval overlap over 10 simulations as a function of system 
size TV, for connectivity degree distribution p{k) ~ k^^ (with p{k) = for A; < 3, at 
T = 1 and p — 5). The dotted line corresponds to the value predicted by our RS 
population dynamics. 

our theory claims that for p{k) ~ k^'^ with 7 < 3 one will have retrieval at any T. For 
large values of p and T, however, the overlaps, although indeed nonzero, become smaller 
and hence our numerical accuracy decreases. Moreover, the equilibration times in the 
simulations grow rapidly as p increases. The discrepancies at high temperatures and 
large p are, we believe, due to finite size effects. In figure 121 we show that the agreement 
between theory and simulations indeed improves for larger system sizes. 

6. Conclusions 

We have solved attractor neural network models on random graphs with arbitrary 
connectivity distributions p{k), using the replica method within RS ansatz, in the spirit 
of [221 1221 HZ]- The RS order parameters are the effective cavity field distributions 
W^{h) in each sublattice, or equivalently, the distributions of messages Q^{u). Second 
order phase transitions from the paramagnetic (P) phase to a retrieval (R) or spin glass 
(SG) phase could be derived explicitly, given the assumption that these transitions 
are second order and provided the second moment {k"^) of the connectivity degree 
distribution of the graph is finite. The overlap and spin glass order parameters in 
each phase can in principle be calculated via a population dynamics algorithm. The 
latter is limited by numerical accuracy when the values of these order parameters are 
small (as for large T and a). We find that the retrieval region in the phase diagram 
is larger for fat-tailed degree distributions than for those with exponential decay (e.g. 
Poisonian), but it is not clear whether this can be exploited in associative memories 
since it goes at the cost of the magnitude of the retrieval overlaps. The possible 
occurrence of replica symmetry breaking is beyond the scope of this paper. Within 
our numerical accuracy, we can conclude that upon comparing the results of our replica 
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symmetric theory (including population dynamics) to numerical simulations, for degree 
distributions ~ [k{k-\-l){k-\-2)]~^ and p{k) ~ /c"''', we obtain satisfactory agreement. 
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